%
% mg_matrix_init.m
%
%   initialize the matrices used in the mulitgrid algorithm.
%
function M = mg_matrix_init(A,levels,tag,nVoffset);
% tag can be 'u','v', or 'p'
% nVoffset gives the flexibility of node-centering or cell-centering.
%  e.g. nVoffset=-1 :  the nodes excluding the boundaries are variables
%  e.g. nVoffset=+1 :  the nodes including the boundaries are variables

  % initialize the data structure
  %
  M.A  = cell(levels,1);
  M.Ir = cell(levels,1);
  M.Ip = cell(levels,1);

  if nargin<4
      nVoffset = 0;
  end
  
  % initialize 
  %
  Ac = A;
 
  for k=levels:-1:1
    N = 2^k+nVoffset;  % size of the current grid
     
    % set the operator
    %
    M.A{k}  = Ac;

    % transfer matricies
    %
    if tag == 'u6'
        interp  = interp_matrix_u_6pt(N,N);
    elseif tag == 'u1'
        interp  = interp_matrix_u(N,N);
    elseif tag == 'v6'
        interp  = interp_matrix_v_6pt(N,N);
    elseif tag == 'v1'
        interp  = interp_matrix_v(N,N);
    elseif tag == 'p0'
        interp  = interp_matrix_p_const(N,N);
    elseif tag == 'p1'
        interp  = interp_matrix_p(N,N,nVoffset);
    else
        error('unsupported variable tag!');
    end

    M.Ip{k} = interp;
    % the denominator 4 follows from the dimensionality being 2
    %  and the coefficients of the bilinear interpolation.
    % In other words, we sum up the coefficients of fine values
    %  that are influenced from a single coarse value.
    M.Ir{k} = interp'/4;
    
    % compute the new coarse grid matrix
    %
    Ac = M.Ir{k} * Ac * M.Ip{k};
  
  end
